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Interpolation  of  the  Radial  Velocity  Data 
From  Coastal  Hf  Radars 
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Abstract 


In  recent  years,  monitoring  nearshore  surface  currents  became  an  important  appli¬ 
cation  of  the  high-frequency  radar  (HFR)  technology.  The  Doppler  shifts  of  backscat- 
tered  radio  signals  from  surface  waves  provide  the  surface  velocity  component  in  the 
direction  of  a  HFR  beam.  These  radial  velocities  observed  by  multiple  (usually  two) 
radars  have  to  be  combined/intcrpolated  to  produce  the  gridded  vector  field,  which  can 
be  used  in  applications.  In  view  of  a  relatively  high  (5-10  cm/s)  HFR  measurement 
errors  of  the  radial  velocities,  interpolation  algorithms  which  take  into  account  addi¬ 
tional  constraints  on  the  velocity  field  (such  as  those  imposed  by  the  coastlines  and 
model  dynamics)  are  of  particular  value. 

In  this  chapter,  recent  developments  in  the  radial  velocity  processing  methods  are 
reviewed.  The  topics  include  advanced  optimal  interpolation  techniques,  kinemati¬ 
cally  constrained  Galcrkin  and  the  2d  variational  interpolation  methods,  and  the  dy¬ 
namically  constrained  assimilation  of  the  HFR  data  using  numerical  models. 

Accurate  monitoring  of  the  velocity  field  may  suffer  from  occasional  malfunction 
of  a  radar  which  causes  a  substantial  data  loss  on  a  relatively  short  (3-30  hours)  time 
interval.  We  examine  performance  of  the  gap-filling  technique  based  on  empirical 
orthogonal  function  analysis  of  the  radial  velocity  observations  and  demonstrate  its 
performance  in  tidally-driven  coastal  environments. 

1.  Introduction 

The  technology  of  monitoring  near-coastal  currents  by  High  Frequency  Radars  (HFRs)  has 
been  rapidly  developing  tn  the  past  decade.  HFR  observations  are  now  extensively  used 
to  study  near-shore  circulation  under  a  large  variety  of  environmental  conditions  (e.g.,  [40, 
16,  17,  9,  43,  45,  1 1])  helping  to  solve  many  applied  problems  in  the  coastal  regions.  At 
present,  the  HFR  surface  current  mapping  technology  is  capable  to  provide  surface  velocity 
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maps  in  an  area  up  to  200  km  offshore,  with  space-time  resolutions  of  0.2-1  hour/1-10  km 
depending  on  the  particular  system  and  practical  application. 

An  important  question  in  dealing  with  HFR  data  is  the  problem  of  retrieving  the  2d 
velocity  vector  maps  from  the  velocity  components  measured  along  the  radar  beams.  The 
most  commonly  used  technique  (e.g.,  [45])  is  based  upon  local  interpolation  (LI)  of  the 
radial  velocity  data.  The  method  is  formulated  as  an  unweighted  least  squares  minimization 
problem  and  takes  into  account  both  measurement  errors  and  the  effect  of  geometric  dilution 
of  precision  [10]  associated  with  the  angle  between  the  HFR  beams  intersecting  in  a  given 
grid  cell.  The  LI  method,  being  a  particular  case  of  the  more  general  optimal  interpolation 
(01)  technique,  has  been  further  improved  by  employing  more  sophisticated  correlation 
functions  [21,  22]. 

Another  algorithm,  which  has  been  under  development  in  recent  years  [26,  20]  is  the 
open-boundary  modal  analysis  (OMA).  This  technique  can  be  viewed  as  an  extension  of 
the  normal  mode  analysis  [13,  32]  which,  similarly  to  the  later  method  of  Park  et  al.  [39], 
employs  decomposition  of  a  2d  vector  field  into  divergent  and  rotational  components.  Apart 
from  capability  to  avoid  explicit  specification  of  the  poorly  known  error  covariance,  the 
OMA  technique  automatically  takes  into  account  the  kinematic  constraints  imposed  on  the 
velocity  field  by  the  coastlines. 

An  alternative  approach,  proposed  recently  by  Yaremchuk  and  Sentchev  [51,  53],  em¬ 
ploys  two-dimensional  variational  (2dVar)  technique  to  map  HFR  radial  velocity  data  on 
the  gridded  vector  field.  In  contrast  to  the  method  of  Park  et  al.,  [39],  this  algorithm  does 
not  require  specification  of  a  background  velocity  field  and  its  covariance  structure  that  may 
be  lacking  in  many  applications.  Similar  to  OMA,  the  2dVar  method  is  non-local  (the  re¬ 
sult  of  interpolation  at  a  given  grid  point  depends  on  all  the  observed  radial  velocities)  and 
kinematically  constrained  (the  interpolated  velocity  field  is  subject  to  constraints  imposed 
by  the  coastlines).  In  addition,  the  2dVar  method  [51]  allows  to  impose  constraints  on  the 
structure  of  the  divergence  and  vorticity  fields,  that  are  important  in  practical  applications, 
such  as  search  and  rescue,  oil  spill  control  etc. 

Accuracy  in  the  reconstruction  of  the  surface  velocity  field  measured  by  HFRs  can  be 
further  improved  by  synthesis  with  other  observations  and  with  dynamical  constraints  pro¬ 
vided  by  the  numerical  models.  One  of  the  earlier  attempts  of  this  kind  was  made  by  Lewis 
et  al.,  [28],  who  employed  HFR  data  to  constrain  wind  stress  forcing  of  a  primitive-equation 
model.  With  subsequent  development  of  the  data  assimilation  (DA)  techniques,  the  HFR 
data  were  widely  used  to  constrain  numerical  models  in  the  framework  of  both  sequential 
[35,  36,  9]  and  4dVar  [19,  24]  DA  schemes.  Although  being  more  general,  the  dynamically 
constrained  methods  of  HFR  data  interpolation  are  significantly  more  computationally  ex¬ 
pensive  and  often  require  a  lot  of  preliminary  tuning  and  quality  control.  For  that  reason, 
operational  centers  currently  use  more  simple  interpolation  methods. 

An  obvious  advantage  of  HFR  observations  is  their  availability  in  real  time  with  nearly 
continuous  temporal  and  spatial  coverage  fully  compatible  with  the  resolution  of  numerical 
models  of  coastal  circulation.  However,  the  back-scattered  HFR  signals  suffer  from  to 
numerous  distortions  of  artificial  and  natural  origin.  As  a  consequence,  estimates  of  the 
along-beam  sea  surface  velocities  extracted  from  the  Doppler  shifts  of  the  signals  become 
unusable,  resulting  in  numerous  gaps  in  spatial  coverage.  These  gaps  may  strongly  degrade 
performance  of  the  interpolation  algorithms  (e.g.,  [20]). 
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A  natural  way  to  fill  these  gaps  is  to  take  into  account  space-time  correlations  between 
the  radial  velocities.  The  aforementioned  DA  algorithms  is  the  most  straightforward  and 
general  approach.  The  underlying  idea  is  to  combine  dynamical  constraints  of  a  model  with 
the  history  of  dense  spatio-temporal  coverage  by  HFRs  to  produce  the  “best”  estimate  of 
the  surface  velocity  at  a  given  time.  This  approach,  however,  has  a  number  of  drawbacks 
hindering  its  implementation  for  real-time  HFR  data  analysis:  Beyond  a  relatively  high 
computational  cost,  the  dynamically  constrained  DA  schemes  have  a  large  number  of  free 
parameters  whose  statistics  is  poorly  known.  The  most  problematic  among  those  are  the 
open  boundary  conditions,  which  are  the  major  contributors  to  slow  convergence  of  the 
HFR  DA  schemes  which  typically  involve  lengthy  open  boundaries. 

Although  the  2d  01  methods  are  computationally  cheaper  than  the  variational  schemes 
involving  dynamical  information,  they  may  perform  poorly  in  the  presence  of  large  gaps 
in  the  data  because  information  on  the  spatial  structure  of  the  velocity  field  within  the  gap 
is  implicitly  drawn  from  the  idealized  covariance  function,  which  looses  accuracy  at  large 
distances.  A  certain  improvement  of  the  covariance  models  can  be  obtained  by  considering 
their  truncated  expansions  in  the  empirical  orthogonal  functions  (EOFs),  a  technique  suc¬ 
cessfully  used  in  Kalman  filtering  (e.g.,  [49])  and  variational  data  assimilation  [14],  [52]. 
The  EOF-based  estimates  of  the  covariances  rely  upon  time  averaging  and,  therefore,  may 
be  successfully  applied  not  only  to  model  output  but  also  to  data  sets  with  nearly  continu¬ 
ous  space-time  coverage  such  as  sea  surface  temperature  (SST)  or  HFR  data.  Beckers  and 
Rixen  [7]  proposed  an  iterative  EOF-based  technique  for  filling  gaps  in  the  gridded  SST 
images,  which  was  successfully  applied  in  the  Adriatic  [1,  2].  Kondrashov  and  Ghil  [23] 
developed  the  method  further  by  including  time  correlations  under  the  assumption  of  statis¬ 
tical  stationarity  of  the  observed  fields.  Yaremchuk  and  Sentchev  [53]  successfully  applied 
this  gap-filling  (GF)  technique  to  the  HFR  data. 

In  this  chapter  we  give  an  overview  of  the  radial  velocity  interpolation  methods  with 
a  focus  on  the  kinematically  constrained  2d  interpolation  schemes  and  the  GF  techniques, 
that  are  likely  to  be  introduced  operationally  in  the  near  future. 

The  chapter  is  organized  as  follows:  In  the  next  section  we  describe  the  LI,  OMA  and 
2dvar  algorithms  and  compare  their  performance  in  the  framework  of  numerical  experi¬ 
ments  with  simulated  and  real  data.  In  Section  3  we  give  an  overview  of  the  accumulated 
experience  in  the  dynamically  constrained  interpolation  of  the  HFR  data  and  their  synthesis 
with  other  observations.  Section  4  describes  applications  of  the  HFR  GF  algorithm  with 
both  simulated  and  real  data  acquired  in  tidally-driven  coastal  environments. 

2.  Dynamically  Unconstrained  Interpolation  Methods 

The  general  approach  discussed  in  this  section  can  be  classified  as  the  two-dimensional  01, 
or  spline  interpolation  [34]:  The  corresponding  least-squares  algorithms  differ  from  each 
other  by  specifying  either  the  covariance  function  or  its  inverse.  In  application  to  HFR 
data  the  01  algorithms  were  implemented  using  empirical  error  covariances  deduced  from 
“normal  modes”  [32],  “open-boundary  modes”  [20]  and  the  data  [21],  [22].  The  spline 
formulation  [5 1  ]  defines  the  inverse  covariance  by  the  roughening  operator  which  penalizes 
the  grid-scale  variability  in  the  2d  velocity,  divergence  and  vorticity  fields. 
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2.1.  Description  of  the  Methodologies 

2.1.1.  Optimal  Interpolation 

In  application  to  interpolation  of  the  radial  velocities  v,  the  traditional  OI  approach  is  formu 
lated  as  the  following  minimization  problem:  find  the  interpolated  velocity  field  ii(jc)  such, 
that  given  a  set  of  n  HFR  velocity  observations  v,-,  /  =  I,...,/i  around  x ,  the  interpolation 
error  variance  e2  is  minimized: 


(e2)  =  (£w/(v,  -«(xi)  ri)2)  -♦  min  (1) 

f*l  1 

Here  angular  brackets  denote  statistical  average,  r,  stand  for  the  radar  beam  directions  in 
the  observation  points  and  w,  are  the  weights  derived  from  the  statistics  of  u(x)  (e.g.,  [30]). 
The  approach  requires  a  statistical  model  for  the  correlation  tensor  (uu  )  of  the  interpolated 
field,  which  may  not  be  readily  available  from  observations.  Therefore,  in  many  practical 
applications,  eq.  ( 1 )  is  simplified  by  selecting  the  grid  which  is  coarse  enough  to  contain  at 
least  two  observation  points  from  different  radars  within  a  given  grid  cell.  By  neglecting  the 
decay  of  correlations  within  the  cell,  the  problem  can  be  reformulated  locally  (e.g.,  [45]): 
Find  the  velocity  vector  u  at  a  given  grid  point  such  that  deviations  of  its  projections  on 
the  beam  directions  rx  at  the  surrounding  observation  points  are  minimized.  In  the  most 
common  case  with  only  two  radars  illuminating  a  grid  point,  the  solution  is 

u  =  (vi cos(X2  -  V2COSOI1)/ sin((X|  -  0C2)  (2) 

v  =  (v2Sinai  —  vi  sin(X2)/ sin(oi|  -  (X2)  (3) 

where  v\  2  are  the  observed  radial  velocities  and  CX12  denote  angles  between  the  coordinate 
axes  and  the  beam  directions  rlf 2  at  the  grid  point.  This  local  interpolation  (unweighted 
least  squares)  method  has  been  widely  used  in  operational  centers  because  of  its  robustness 
and  simplicity.  However,  due  to  discontinuous  correlation  functions  it  often  yields  spurious 
values  along  baselines  between  HFR  stations.  Moreover,  as  it  is  seen  from  (2-3),  the  algo¬ 
rithm  looses  accuracy  with  the  reduction  of  the  angle  ot|  -  a  2  between  the  radar  beams  in 

the  vicinity  of  the  interpolated  point,  and  it  would  not  work  at  all  with  only  one  radar,  or  in 

the  regions  of  sparse  HFR  coverage. 

2.1.2.  Normal  Mode  Analysis 

Limitations  imposed  by  the  local  nature  of  the  LI  algorithm  can  be  formally  overcome  by 
the  Galerkin  methods  which  employ  global  basis  functions  covering  the  entire  interpolation 
domain.  The  early  technique  of  this  type  was  proposed  by  Eremeev  et  al.  [13]  for  interpo¬ 
lation  of  the  sparsely  sampled  velocities.  Later,  Lipphardt  et  al  [32]  applied  the  technique 
to  HFR  observations.  To  perform  the  interpolation,  the  2d  velocity  field  is  decomposed  into 
divergent  and  rotational  components  defined  by  the  velocity  potential  <p  and  the  stream  func¬ 
tion  \| /.  The  stream  function  and  velocity  potential  are  then  expanded  in  a  predetermined  set 
of  2d  functions  (normal  modes)  with  varying  smoothness  (eigenfunctions  of  the  Laplacian 
operator  A  in  the  interpolation  domain  Q).  Expansion  coefficients  are  then  optimized  to  fit 
the  data. 
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The  normal  mode  method  was  further  developed  by  Lekien  and  Coulliette  [26],  Kaplan 
and  Lekien  [20]  who  supplemented  the  Laplacian  eigenfunctions  by  the  additional  expan¬ 
sion  functions  (open  boundary  modes).  These  functions  were  introduced  to  bring  more  re¬ 
alism  near  the  open  boundary,  and  were  defined  as  solutions  to  the  Poisson  equation  forced 
by  the  Dirichlet  boundary  conditions  of  variable  smoothness.  The  algorithm  is  controlled 
by  two  free  parameters:  the  spatial  length  scale  L  which  defines  the  number  of  modes  used 
for  interpolation,  and  the  regularization  parameter  K  which  penalizes  the  amplitude  of  the 
modes. 

An  advantage  of  the  OMA  method  is  its  ability  to  control  spatial  roughness  by  pre¬ 
scribing  a  limited  number  of  the  smoothest  eigenfunctions  for  interpolation.  Additionally, 
the  method  is  capable  to  constrain  the  interpolation  problem  kinematically  by  taking  into 
account  the  conditions  imposed  by  the  coastlines.  However,  the  OMA  method  lacks  flex¬ 
ibility  in  representing  localized  small-scale  features  as  well  as  flow  structures  near  open 
boundaries,  where  the  rotational  component  of  the  velocity  field  is  assumed  to  be  tangent 
to  the  boundary.  Besides,  certain  difficulties  may  arise  when  dealing  with  gappy  data,  es¬ 
pecially  when  the  horizontal  size  of  a  gap  is  larger  than  the  minimal  length  scale  resolved 
[20].  These  drawbacks  of  the  OMA  technique  may  impose  certain  limitations  on  taking  the 
full  advantage  of  the  HFR  data  whose  spatial  resolution  strongly  varies  across  the  domain. 

2.L3.  Variational  Method 

Interpolation  of  the  HFR  data  with  variational  technique  is  more  flexible,  as  it  is  not  con¬ 
strained  by  a  predetermined  set  of  functions.  Instead,  the  interpolating  functions  in  use 
(eigenvectors  of  the  Hessian  matrix)  depend  both  on  the  shape  of  the  coastline  and  spatial 
distribution  of  the  data  points.  The  latter  feature  (lacking  in  OMA)  is  especially  useful  in 
practice,  because  HFR  data  often  have  gaps  in  spatial  coverage,  and  the  issue  of  filling  those 
gaps  is  important. 

Another  useful  property  of  the  2dVar  method  is  that,  similar  to  LI,  it  operates  in  the 
velocity  space,  and  thus  requires  less  sophisticated  operators  for  projection  of  the  unknown 
gndded  velocities  on  the  radial  components  of  the  flow  speed.  This  property  provides  better 
conditioning  of  the  Hessian  matrix,  faster  convergence  of  the  minimization  algorithm  and 
improved  computational  efficiency. 

The  basic  formulation  of  the  2dVar  scheme  [51,  53]  is  the  following:  Find  the  velocity 
field  «(jc)  such,  that  the  cost  function 

7==  2^  E  v*]2  +  ^/ [w</(Ad,v“)2'f  VVC(Acur,w)2  +  'v"(Aw)2]  —  m;0  (4) 

*-i  n 

is  minimized  with  respect  to  the  grid  point  values  u  under  the  constraint  //(dD)  =  0.  Here 
K  is  the  number  of  HFR  observations  in  f2,  A  =  / dQ  and  is  the  local  interpolation 
operator  which  projects  the  unknown  velocity  vectors  onto  the  kih  observation  point  from 
the  apexes  of  the  grid  cell,  enveloping  that  point.  Factors  C^2,Wd,Wc  and  Wu  are  the 
inverse  error  variances  of  the  corresponding  squared  quantities,  so  that  J  could  be  treated  as 
the  argument  of  the  Gaussian  pdf  T(u )  defined  on  the  2Af*dimensional  space  of  the  gridded 
velocity  fields  u:  t P(w )  ~  exp[-7]. 

As  it  is  seen  from  (4),  regularization  of  the  interpolation  problem  is  performed  by  pe¬ 
nalizing  not  only  the  grid-scale  components  of  u  (the  last  term),  but  also  of  its  curl  and 
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Figure  1.  Setting  of  the  numerical  experiments:  Left  panel  shows  the  ’’true”  velocity  field 
u*  which  has  nine  eddy  structures  superimposed  on  two  jets:  one  following  the  coastline 
and  another  located  farther  offshore.  The  right  panel  shows  the  vorticity  field  curia7.  Three 
coastal  radars  sample  the  radial  velocities  along  beam  directions  binned  at  2  km  radial 
and  5°  azimuthal  resolution.  Sampling  points  are  shown  by  gray  dots.  Gray  rectangle  in 
the  upper  right  panel  envelops  simulated  gap  in  the  HFR  data.  Radar  positions  are  shown 
by  black  dots  on  the  coastline.  The  interpolation  domain  and  the  coastline  are  similar  to 
the  ones  used  in  [20]  for  OMA  processing  of  HFR  observations  in  the  Bodega  Bay.  The 
coordinate  system  is  rotated  clockwise  (north  is  on  the  right).  Contour  interval  for  vorticity 
is  5*  10"5  s"1. 

divergence  (terms  weighted  by  Wd  and  Wc).  This  is  done  in  an  attempt  to  retrieve  the 
larger-scale  components  of  divw  and  curlw  that  are  important  in  applications. 

2.2.  Comparison  of  the  Interpolation  Techniques 

Performance  of  the  interpolation  methods  outlined  in  the  previous  section  are  compared 
using  two  complimentary  techniques.  First,  the  HFR  data  are  simulated  using  a  realistic 
geometry  of  the  3-radar  experiment  in  the  Bodega  Bay  conducted  in  spring  and  summer  of 
2003  (Fig.  1).  In  the  second  series  real  data  acquired  on  July  30,  2003  in  the  Bodega  Bay 
and  offshore  the  coast  of  Brittany  has  been  used. 

2.2.1.  Experiments  with  Simulated  Data 

Simulated  data  experiments  are  performed  by  varying  the  interpolation  method  (2dVar,  LI 
and  OMA)  and  coverage  of  the  domain  by  observations  (with  and  without  the  gap  shown 
in  the  right  panel  of  Figure  1).  Within  each  series  of  experiments  the  noise  level  v  in  the 
simulated  observations  was  also  varied:  the  radial  velocities  v*k  observed  at  points  jr*  were 
defined  by  adding  white  noise  w  to  projections  of  the  true  currents  ul  on  the  beam  directions 

rk- 

vl  =  {PkJu'ij-rk)+Ww.  (5) 

Here  V  is  the  typical  magnitude  of  ul  and  v  is  the  scalar  parameter  whose  reciprocal  has 
the  meaning  of  signal/noise  ratio.  Three  values  of  v  (0.1,  0.3  and  0.5)  were  tested  within 
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each  series.  In  correspondence  with  the  equation  (4),  HFR  measurement  errors  in  (5)  were 
defined  as  o(vj)  =  W. 

To  assess  the  impact  of  gaps  in  the  spatial  coverage  by  HFR  observations,  two  simulated 
data  sets  were  used:  One  with  the  full  coverage  by  the  radial  velocity  data  and  another 
without  observations  in  the  rectangular  domain  (right  panel  in  Fig.l).  These  data  sets  had 
K- 201 1  and  K- 1699  observation  points  respectively. 

The  quality  of  interpolation  was  monitored  by  several  parameters.  Misfit  between  the 
interpolated  field  and  the  data  was  defined  as  mi j  —  |v*  -  vj|/|vj|,  where  v *  is  the  projec¬ 
tion  of  interpolated  velocity  on  the  radar  beam  at  the  kth  measurement  point  and  overline 
denotes  averaging  over  the  data  points.  Velocity  error  ev  was  defined  as  the  mean  absolute 
difference  between  the  true  u!  and  interpolated  u  currents  normalized  by  V: 

ev=  (\u'-u\)/V,  (6) 

where  angular  brackets  denote  averaging  over  the  interpolation  grid.  Similar  expressions 
are  used  to  assess  the  interpolation  qualities  e<j,  ec  of  the  divergence  and  vorticity  fields: 

erf<|div(i/'  -i/)|)/(|div7/'|);  ec<|curl(i/'  -  w)|)/<|curlw'|)  (7) 

To  simplify  the  analysis,  the  inverse  vanance  Wu  is  set  to  zero  in  all  the  experiments.  It 
is  also  assumed  that  Wd  =  y  2Wcy  where  y  =  0.2  has  the  meaning  of  the  Rossby  number 
in  most  situations,  the  2d  vorticity  field  is  much  stronger  than  divergence,  as  the  latter  is 
driven  by  relatively  weak  vertical  motions  acting  against  gravity.  Finally,  since  the  result  of 
interpolation  depends  only  on  the  relative  magnitudes  of  the  terms  in  the  cost  function,  it 
is  convenient  to  charactenze  the  2dVar  algorithm  by  the  single  non-dimensional  parameter 
Wf  =  v2V2WcK/4A8xa  =  \2KLA /AAhx2  which  has  the  meaning  of  the  ratio  between  the 
vorticity  regularization  and  the  data  misfit  terms.  Here  6x  =  2  km  is  the  grid  step  and  L  =  5.9 
km  is  the  typical  horizontal  scale  of  the  currents  shown  in  Fig.  1.  Since  Wd  =  y~  Wc, 
the  value  of  Wf  entirely  defines  the  shape  of  the  cost  function  for  a  given  field  and  data 
configuration.  More  details  on  the  experimental  setting  can  be  found  in  [5 1 J. 

Results  of  the  experiments  with  three  radars  are  illustrated  in  Table  1  and  Figures  2  and 
3.  Table  1  contains  errors  in  approximation  of  the  true  field  and  the  data  for  different  noise 
levels.  It  is  seen  that  the  2dVar  algorithm  provides  better  quality  of  reconstruction  than 
LI  and  OMA  methods,  although  the  reconstructed  velocity  fields  look  visually  similar  (left 
panels  in  Fig.  2).  In  that  sense  the  capability  of  HFR  measurements  to  capture  the  structure 
of  the  flows  in  Fig.  1  is  quite  remarkable.  With  the  increase  of  observation  noise  all  the 
algorithms  loose  precision,  but  still  provide  a  noise-consistent  fit  to  the  data  in  the  entire 
range  of  v  (cf.  columns  2-3  and  6-7  in  Table  1).  The  only  exception  is  the  OMA  algorithm 
forv=0.1  (second  line  in  Table  1). 

Reconstruction  of  the  divergence  and  vorticity  fields  is  less  accurate  since  horizontal 
derivatives  amplify  grid-scale  noise.  This  is  especially  visible  in  the  behavior  of  because 
the  divergence  field  was  set  to  be  five  times  weaker  than  vorticity  by  the  design  of  the 
experiments.  Nevertheless,  in  contrast  to  the  OMA  and  LI  methods,  the  2dVar  algorithm 
captures  signs  and  positions  of  the  major  structures  in  the  divergence  field  at  the  “practical” 
noise  levels  of  0.1  and  0.3  (upper  right  panel  in  Figures  2). 

Compared  to  the  divergence,  the  vorticity  field  is  reconstructed  with  much  higher  qual¬ 
ity:  Both  the  amplitudes  and  locations  of  most  of  the  structures  are  well  reproduced  even 
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Table  1.  Dependence  of  the  interpolated  field  error  parameters  md,  ev,  ec,  and  e<j  on  the 
reconstructed  velocity  field  and  noise  level  v  in  the  simulated  HFR  data.  The  2dVar, 
OMA  and  LI  results  are  shown  respectively  in  the  upper,  middle  and  lower  lines  of 
the  table  cells.  The  smallest  errors  are  boldfaced.  Note  that  LI  errors  were  computed 
over  subdomains  which  did  not  include  near-coastal  regions  and  the  gap  (see  Fig.  2- 
3).  Parameters  for  the  2dVar  and  OMA  interpolation  schemes  are  shown  in  the  right 
column.  The  2dVar  experiments  with  v=0.5  were  made  under  the  divergence-free 
approximation  (W?  =  106) 
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0.56 
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6,  10  3 

0.34 

0.41 

0.57 

2.99 

0.49 

0.50 

0.69 

2.36 

at  the  noise  level  of  v=0.3  (cf.  right  panel  in  Fig.l  and  left  panels  in  Fig.  3  ).  The  gap  in 
data  coverage  is  also  handled  well:  the  overall  increase  in  of  the  interpolation  errors  eVy  e<j 
and  ec  in  Table  1  is  consistent  with  the  fraction  of  the  interpolation  domain  occupied  by 
the  gap,  whereas  the  saddle-like  structure  of  the  vorticity  field  in  the  gap  is  reconstructed 
quantitatively  by  the  2dVar  method. 

The  velocity  patterns  generated  by  the  three  methods  (left  panels  in  Figure  2)  differ 
in  terns  of  the  velocity  interpolation  error  ev.  The  difference  becomes  more  evident  after 
taking  the  divergence  of  the  fields  (right  panels  in  Figures  2  and  3):  Since  OMA  does  not 
impose  any  smoothness  constraint  on  divu,  the  corresponding  divergence  field,  although 
being  3  times  weaker  than  vorticity,  looks  rather  chaotic  with  the  formal  error  e^=1.52.  A 
similar  featureless  pattern  is  produced  by  LI,  with  the  difference  that  divergence  in  near¬ 
coastal  areas  cannot  be  estimated  at  all  due  to  either  single-radar  measurements  or  to  beam¬ 
crossing  angle  limitation  (nearly  parallel  radar  beams). 

Introduction  of  the  gap  in  data  coverage  enhances  the  difference  between  the  three  meth 
ods  (cf.  columns  4-8  and  5-9  in  Table  1  keeping  in  mind  that  LI  errors  are  computed  outside 
the  gap).  A  probable  reason  for  the  difference  between  OMA  and  2dVar  is  emergence  of 
the  spurious  maxima  in  both  vorticity  and  divergence  fields  inside  the  gap  in  the  OMA  case 
(cf.  upper  and  middle  panels  in  Fig.  3).  In  the  2dVar  formulation,  vorticity  and  divergence 
fields  cannot  have  local  maxima  inside  data-void  regions,  because  their  variation  within  the 
gap  closely  approaches  the  behavior  of  a  harmonic  function. 

Comparison  of  the  2dVar,  OMA  and  LI  lines  in  Table  1  shows  that  OMA  code  tends  to 
provide  less  precise  fit  to  the  data  at  10%  noise  level  without  the  gap  (line  2  in  Table  1). 
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Figure  2.  Comparison  between  the  three  methods  of  interpolation.  The  noise  level  v  and 
interpolation  errors  are  shown  in  the  upper  right  comers  of  the  corresponding  panels. 
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re  3.  Comparison  between  the  interpolations  of  three  radars  with  a  gap  in  observa- 
.  The  noise  level  v  and  interpolation  errors  are  shown  in  the  upper  right  comers  of  the 
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This  can  be  partly  explained  by  a  relatively  low  number  of  degrees  of  freedom  (number  of 
modes)  involved.  In  the  OMA  formulation  the  number  of  modes  is  proportional  to  ( D/L )2 
[20],  where  D  is  the  horizontal  size  of  the  domain  and  L  is  the  cut-off  length  scale.  For  the 
reported  experiments,  however,  the  optimal  length  scale  L  was  close  to  6  km  (right  column 
in  Table  1 ).  Such  length  scale  corresponds  to  200-230  eigenmodes,  whose  amplitudes  were 
varied  to  fit  the  data  in  201 1  (1699  with  the  gap)  observation  points.  In  an  attempt  to  achieve 
a  better  fit  at  low  noise  levels,  we  tried  to  increase  the  number  of  modes  by  reducing  L  to 
2-3  km,  but  that  required  an  increase  of  regularization  parameter  K,  otherwise  interpolation 
patterns  appeared  too  noisy,  possibly  because  of  the  ill-conditioning  of  the  system  matrix. 

Overall,  Table  1  shows  that  the  2dVar  method  performs  similar  to  LI  and  better  than 
OMA  at  v=0. 1  and  somewhat  better  at  v=0. 3-0.5.  When  the  gap  is  present  in  the  data  (lines 
1-2  and  4-5  in  Table  1)  2dVar  keeps  a  significant  advantage  to  OMA  up  to  v=0.3.  Both 
non-local  methods  (OMA  and  2dVar)  are  better  than  LI  because  of  their  ability  to  estimate 
currents  within  the  gap  and  close  to  the  coastline. 


Table  2.  Same  as  in  Table  1,  but  for  the  experiments  with  2-radar  configurations 
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In  practice,  there  are  often  situations  when  a  radar  stops  operating  due  to  hardware 
failure  or  some  other  reasons.  In  such  case  local  interpolation  methods  often  fail  in  a  large 
number  of  gridpoints,  because  they  require  at  least  two  velocity  components  for  retrieving 
the  velocity  vector  in  a  grid  cell.  The  OMA  and  2dVar  algorithms  are  essentially  non-local 
and  therefore  have  an  ability  to  interpolate  the  velocity  field  over  distances  exceeding  the 
grid  cell  size  Sir. 

To  investigate  the  performance  of  the  schemes  in  such  situations,  we  switched  off  the 
rightmost  (northern)  and/or  middle  radars  and  examined  the  interpolation  patterns  both  with 
and  without  the  gap  in  the  data.  These  experiments  also  allowed  us  to  assess  the  accuracy 
of  interpolation  in  the  regions  where  data  density  was  less  or  close  to  one  observation  per 
grid  cell:  After  removing  the  northern  radar  such  regions  emerge  in  the  upper  (western)  and 
right  (northern)  parts  of  the  domain. 
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Figure  4.  Comparison  between  the  interpolations  with  two  radars.  The  noise  level  v  and 
interpolation  errors  are  shown  in  the  upper  right  comers  of  the  corresponding  panels. 

Figure  4  gives  an  indication  that  OMA  algorithm  is  less  accurate  than  2dVar  in  such 
regions:  A  visual  comparison  of  the  vorticity  distributions  with  those  in  Figure  1  shows  that 
OMA  errors  tend  to  be  larger  when  y  >  40  km  or  x  >60  km.  The  LI  algorithm  performs 
much  worse:  derivatives  of  the  velocity  field  cannot  be  estimated  not  only  in  the  region  of 
single-radar  coverage  (x  >60  km)  but  also  near  the  coastline.  Quantitative  assessment  of 
ec  for  the  three  interpolation  patterns  in  Figure  4  also  shows  that  the  regions  of  low  data 
density  are  the  major  contributors  to  the  larger  value  of  ec  for  the  OMA  pattern  (0.45  vs 
0.34  in  Figure  4).  The  value  of  ec  for  LI  (0.39)  cannot  be  objectively  compared  with  these 
numbers  because  it  was  computed  by  averaging  over  much  smaller  area. 

Figure  5  shows  an  example  of  interpolation  with  two  radars  and  the  gap  in  the  data.  The 
difference  between  2dVar  and  OMA  is  already  evident  from  the  velocity  patterns:  OMA 
produces  a  spurious  jet  within  the  gap  which  destroys  two  eddies  at  the  upper  and  right 
edges  of  the  data-void  region.  2dVar  preserves  these  eddies  and  the  saddle-like  structure 
of  the  currents  within  the  gap.  Comparison  of  the  vorticity  fields  (right  panels  in  Figure  5) 
shows  that  OMA  again  produces  a  maximum  inside  the  gap.  The  2dVar  pattern  appears  to 
be  unrealistic  inside  the  gap  as  well,  but  compared  to  OMA  has  more  reasonable  structure 
near  the  gap's  boundary. 

The  overall  results  (Table  2)  indicate  that  in  the  case  of  two-radar  configurations  non¬ 
local  methods  (OMA  and  2dVar)  have  a  considerable  advantage  over  LI  with  2d  Var  showing 
somewhat  better  performance  than  OMA. 

2.2.2,  Real  Data  Examples 

To  assess  the  performance  of  the  interpolation  with  real  data,  velocity  field  from  real  HFR 
observations  in  the  Bodega  Bay  have  been  reconstructed  (Fig.  6).  The  interpolation  pa¬ 
rameters  for  OMA  were  L= 5  km  and  K=I0“4  as  in  [20].  The  2d  Var  algorithm  was  used  in 
divergence-free  approximation  (W/  =  106),  because  at  the  estimated  noise  level  of  v=0.35, 
velocity  scale  of  L=4.4  km  and  sampling  discretization  of  2  km,  it  is  hard  to  obtain  statisti¬ 
cally  confident  estimates  of  the  divergence.  With  1401  observation  points  in  use,  W*c=0.5. 

The  2dVar  and  OMA-  generated  patterns  are  qualitatively  similar,  although  OMA  ve¬ 
locity  is  more  smooth  and  characterized  by  somewhat  larger  misfit  with  the  data.  The  ma¬ 
jor  difference  between  the  2dVar-  and  OMA-generated  patterns  is  observed  in  two  regions 
shown  in  Fig.  6  by  gray  rectangles.  In  the  bay  between  the  southern  (left)  and  middle  radars 
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.  Comparison  between  the  interpolations  of  two  radars  with  a  gap.  The  noise  level 
terpolation  errors  are  shown  in  the  upper  right  comers  of  the  corresponding  panels. 
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Figure  6.  Velocity  field  in  the  Bodega  Bay  on  July,  30,  2003  obtained  by  two  interpolation 
methods.  Misfit  m<j  with  the  radial  velocity  data  is  shown  in  the  upper  right  comers. 


the  2dVar  pattern  shows  an  indication  of  anticyclonic  circulation,  whereas  OMA  produces 
a  broad  offshore  current  there.  The  values  of  ev  computed  by  averaging  over  265  observa¬ 
tion  points  in  that  region  are  0.26  for  2d Var  and  0.27  for  OMA  respectively.  The  region  on 
the  right  should  be  considered  as  a  gap,  since  it  contains  only  4  data  points  near  its  right 
boundary. 

The  LI  interpolation  (not  shown)  appeared  to  have  much  lower  quality  compared  to  the 
results  of  non-local  interpolation.  This  was  evident  not  only  in  terms  of  the  larger  tfv=0.29, 
but  also  visually:  velocity  estimates  do  not  exist  in  the  sparsely  covered  right  comer  of  the 
domain  and  are  too  noisy  near  the  coast  where  radar  beams  are  close  to  parallel. 
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Figure  7.  Surface  currents  in  the  Iroise  Sea  (Northern  Bay  of  Biscay)  on  September  5,  2005 
(02:48GMT)  obtained  by  the  2d  Var  (left)  and  OMA  methods.  50  and  100m  isobaths  are 
shown  by  gray  contours.  HFR  location  at  Point  Garchine  is  shown  by  the  black  dot. 

Figure  7  shows  a  comparison  of  the  surface  current  velocities  obtained  by  the  OMA 
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and  2dVar  methods  on  September  5,  2005  in  the  Iroise  Sea  (northern  Bay  of  Biscay).  This 
area  is  continuously  monitored  by  two  high-frequency  Wellen  Radars  (WERA)  operating 
at  12.4  MHz  since  July  2006.  The  radar  sites  are  located  at  the  Point  Garchine  (Fig.  7) 
and  Brezellec  (not  shown,  5  km  east  of  the  SE  comer  of  the  domain).  The  Beam  forming 
(BF)  method,  actually  used  for  processing  the  radar  data  in  near  real-time,  provides  radial 
velocities  of  surface  current  along  beams  with  3dB  width  of  9°.  The  along-beam  resolution 
is  1.5  km  and  the  maximum  range  is  of  the  order  of  140  km  [33].  Radial  velocities  from 
two  radars  are  interpolated  routinely  with  LI  algorithm  to  provide  surface  current  maps  at 
20-mi n  acquisition  rate  (http://www.previmer.org/observations/courants/radarJifJroise). 

Regional  currents  are  quite  strong  (up  to  4  m/s),  making  it  difficult  to  observe  them  by 
bottom-mounted  currents  meters.  The  flow  structure  is  characterized  by  the  strong  horizon¬ 
tal  shear  (up  to  10  3  s_l)  in  the  vicinity  of  the  Ushant  Islands  and  significant  divergence 
induced  by  small-scale  topographic  features,  controlling  the  velocity  field. 

Comparison  of  the  velocity  patterns  in  Fig.  7  shows  better  flexibility  of  the  2dVar  al¬ 
gorithm  in  fitting  the  data  when  the  local  scales  of  variability  abruptly  change  within  the 
interpolation  domain  as  it  is  seen  on  the  transition  from  the  southwestern  part  of  the  re¬ 
gion  to  the  region  of  Ushant  and  Molene  Islands  in  the  northeast.  The  2dVar  algorithm 
provides  a  significantly  better  resolution  in  that  area,  generating  a  strongly  sheared  flow 
around  the  Ushant  Island  and  significantly  smaller  interpolation  errors  in  this  region:  the 
root-mean-square  deviations  of  the  interpolated  velocities  from  the  observed  radial  veloci¬ 
ties  are  respectively  0.24  and  0.39  m/s  for  the  2dVar  and  OMA-generated  patterns.  Specific 
features  of  regional  circulation,  such  as  the  control  of  the  flow  by  bathymetry,  current  in¬ 
tensification  and  deflection  of  the  currents  north  of  Ushant  Island,  are  better  reproduced  by 
the  2dVar  algorithm  [44]. 


3.  Dynamically  Constrained  Interpolation  Techniques 

Interpolation  techniques  discussed  in  the  previous  section  are  essentially  two-dimensional, 
i.e.  interpolated  velocity  patterns  at  different  times  are  not  connected  neither  statistically  nor 
dynamically.  In  the  last  decades  advanced  interpolation  methods  have  been  developed  that 
produce  dynamically  and  statistically  consistent  estimates  of  the  entire  ocean  state.  These 
methods  combine  dynamical  constraints  from  numerical  models  of  the  oceanic  circulation 
with  statistical  information  from  observations.  By  using  such  approach  it  becomes  possible 
not  only  to  interpolate  the  observed  data  in  space  and  time,  but  also  infer  information  on 
unobserved  variables,  such  as  wind  stress,  bottom  friction  coefficient,  vertical  velocity  etc. 
Another  advantage  of  this  approach  is  its  ability  of  natural  blending  HFR  data  with  other 
types  of  observations.  Combining  numerical  models  and  data  to  make  an  optimal  analysis  is 
referred  to  as  data  assimilation  (DA)  which  is  now  widely  used  in  oceanography.  The  data 
assimilation  methods  are  divided  into  two  categories:  sequential  and  fully  four-dimensional 
in  space-time.  In  this  section,  we  give  a  brief  overview  of  the  existing  applications  of  the 
DA  techniques  to  HFR  data. 
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3.1.  Sequential  Methods 

The  most  widely  used  DA  schemes  optimize  model  solutions  sequentially  in  time:  After 
a  model  forecast  is  adjusted  to  observations,  the  analysis  (adjusted  state)  is  propagated 
forward  in  time  by  the  model  until  new  data  are  available.  The  adjustment  (interpolation)  is 
usually  formulated  as  a  weighted  least  squares  problem  with  the  weights  defined  by  the  error 
covariance  matrix  of  the  observations.  In  the  simplest  DA  scheme,  the  error  covariance  does 
not  change  in  time,  and  the  problem  is  reduced  to  a  sequence  of  OIs  and  model  integrations 
of  the  interpolated  states.  The  approach  when  both  the  state  and  its  error  covariance  are 
propagated  by  the  model  is  referred  to  as  Kalman  filtering  (KF). 

Although  HFR  measurements  represent  a  very  valuable  data  set  for  coastal  ocean  state 
estimation,  DA  of  the  HFR  observations  is  only  at  the  very  beginning  of  the  massive  oper¬ 
ational  use.  Comparatively  simple  techniques,  mainly  based  on  nudging  (e.g.  [50],  [28]) 
allowed  Lewis  et  al.  [28]  to  optimize  poorly  known  atmospheric  forcing  in  the  Monterey 
Bay,  California.  An  approach  is  taken  in  which  the  HFR  data  act  as  if  there  were  an  ad¬ 
ditional  layer  of  water  overlying  the  ocean  surface.  A  pseudo-shearing  stress  resulting 
from  the  difference  between  the  model-predicted  velocity  and  the  Doppler  radar  velocity  is 
added  to  that  of  the  wind  in  order  to  force  a  model.  However,  the  method  had  not  used  any 
knowledge  of  the  error  covariances  that  could  improve  blending  of  the  model -predicted  and 
observed  surface  currents.  Furthermore,  wind-forcing  was  corrected  only  at  model  loca¬ 
tions  that  were  within  the  HFR  footprint,  thus  introducing  a  risk  of  unrealistic  irregularities 
at  the  edges  of  the  footprints. 

Using  a  similar  approach,  Wilkin  et  al.  [50]  performed  a  vertical  extrapolation  of  the 
HFR  surface  velocities.  In  the  experiments,  a  statistical  projection  scheme  based  on  corre¬ 
lations  between  the  HFR  and  the  moored  Acoustic  Doppler  Current  Profiler  data  was  used. 
The  experiments  showed  that  HDR  data  provides  an  invaluable  source  of  information  for 
coastal  ocean  DA  forecast  systems  and  that  even  a  simple  statistically  based  vertical  extrap¬ 
olation  significantly  improves  the  skill  in  estimating  subsurface  velocities  in  New  Jerseys 
shelf  region. 

In  a  number  of  recent  studies  [9,  36,  37,  5]  sequential  DA  schemes  were  used  assimilate 
the  surface  currents  with  covariance  estimates  based  on  the  ensembles.  Oke  et  al.  [36] 
performed  an  assimilation  of  low-pass  filtered  surface  velocity  observations  from  coastal 
HFR  arrays  into  a  primitive  equation  coastal  ocean  model.  The  error  covariance  model  was 
obtained  from  the  ensemble  of  model  simulations  and  involved  correlations  between  the 
surface  velocity  and  the  subsurface  fields  of  temperature,  salinity  and  velocity.  Using  these 
inhomogeneous  error  covariances  among  the  model  variables,  the  3d  model  corrections 
based  on  the  surface-only  HFR  data  were  derived  and  successfully  used  in  the  DA  scheme. 
The  correlation  between  subsurface  current  measurements  and  subsurface  currents  obtained 
from  model  was  improved  two  times  (from  0.42  to  0.78),  when  radar  data  were  assimilated, 
thus  demonstrating  the  effectiveness  of  such  DA  approach. 

In  a  similar  manner,  a  3d  covariance  matrix  for  the  Princeton  Ocean  model  model  has 
been  constructed  by  Breivik  and  Saetra  [9]  who  used  a  statistically  representative  reference 
model  run.  The  matrix  has  been  kept  fixed  throughout  the  assimilation  period.  The  ensem¬ 
ble  KF  based  optimal  interpolation  of  HFR  velocities  recorded  off  the  Norway  coast,  has 
been  used  to  provide  dynamically  balanced  current  fields  generated  by  a  system  of  nested 
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ocean  models  in  near-real  time.  The  system  was  found  to  yield  good  analyses,  whereas  the 
short  range  forecasts  were  significantly  improved  by  assimilating  HFR  data. 

Another  sequential  DA  scheme  based  on  3dVar  approach  was  recently  developed  by  Li 
et  al.  [29]  for  operational  assimilation  of  multiple  data  sets  collected  in  Southern  California 
into  the  Regional  Ocean  Modeling  System  (ROMS).  The  model  is  able  to  assimilate  surface 
current  observations  from  HF  radar  network.  Inhomogeneous  and  anisotropic  error  covari¬ 
ances  have  been  constructed  by  blending  ensemble  runs  with  heuristic  covariance  functions, 
and  a  variety  of  observational  data  collected  during  the  Autonomous  Ocean  Sampling  Net¬ 
work  (AOSN)  experiment  (Monterey  Bay,  August  2003)  was  assimilated.  The  model  ve¬ 
locities  were  compared  with  HFR  and  mooring  observations,  demonstrating  a  reasonable 
accuracy  and  indicating  that  the  system  is  capable  of  reproducing  complex  flows  associated 
with  upwelling  and  relaxation,  as  well  as  the  rapid  transitions  between  them. 

The  approach  of  Lewis  et  al.  [28]  was  recently  revisited  by  Barth  et  al  [4,  5, 6]  who  as¬ 
sumed  that  the  larger  part  of  the  model  error  in  surface  currents  might  be  attributed  to  errors 
in  surface  winds.  Using  the  model  error  covariance  estimated  from  an  ensemble  of  model 
runs  driven  by  different  wind-forcing,  Barthe  et  al.  [5]  assimilated  radial  velocities  recorded 
by  the  HFRs  on  Florida  shelf.  The  assimilation  skill  was  assessed  by  comparing  the  model 
results  with  independent  ADCP  measurements.  As  expected,  the  largest  improvement  was 
observed  at  the  surface  but  the  model  skill  relative  to  the  free  model  run  was  significantly 
improved  also  at  depth.  Later  on,  sequential  assimilation  of  HFR  data  performed  by  Barth 
et  al.,  [6]  improved  the  wind  forcing  instead  of  directly  modifying  the  model  state  vector. 
The  authors  used  an  ensemble-based  assimilation  scheme  and  3d  coastal  circulation  model 
applied  this  time  to  the  German  Bight  region. 

Sequential  methods  described  above  use  a  prior  estimate  of  the  forecast  error  covari¬ 
ance  matrix  that  can  be  difficult  to  obtain  in  practice.  Although  these  methods  produce 
dynamically  and  statistically  consistent  estimates  of  the  ocean  state  and  they  are  capable  of 
improving  the  model  forecast  skills,  they  do  not  produce  a  completely  dynamically  consis¬ 
tent  interpolation  of  the  data  in  both  space  and  time.  This  task  can  be  accomplished  by  the 
4dVar  analysis  and  Kalman  smoothing. 

3.2.  4-Dimensional  Analyses 

Devenon  [12]  was  among  the  first  who  attempted  to  perform  the  dynamically  constrained 
interpolation  of  HFR  data  using  a  variational  approach  and  an  adjoint-based  assimilation 
scheme.  A  simple  numerical  model  (two-dimensional  linearized  spectral  tidal  model)  was 
used  to  smooth  and  interpolate  HFR  observations  in  a  limited  size  coastal  area  (the  Seine 
river  Bay).  The  method  constrains  circulation  to  simultaneously  agree  with  observations 
and  with  the  hydrodynamic  laws  governing  tidal  circulation  in  the  Bay.  Moreover,  the  HFR 
data  have  been  used  to  identify  boundary  conditions  for  the  principal  M2  tidal  constituent, 
and  the  bottom  friction  coefficient. 

Similar  approach  was  used  by  Sentchev  and  Yaremchuk  in  the  Strait  of  Dover  [41] 
and  English  Channel  [43].  The  open  boundary  conditions  of  the  finite-element  spectral 
tidal  model  [27]  were  optimized  to  fit  the  HFR  surface  velocities  and  coastal  tidal  gauge 
data.  The  result  of  interpolation  was  used  for  mapping  of  the  residual  transport  through  the 
Channel,  tidal  dissipation,  and  for  estimation  of  the  energy  flux.  Analysis  of  the  residual 
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flow  stream  function  revealed  a  number  of  permanent  eddies  associated  with  peculiarities  of 
the  coastline  shape  and  inhomogeneities  of  the  bottom  topography.  Extensive  error  analysis 
of  the  results  was  conducted  via  an  explicit  inversion  of  the  Hessian  matrix,  associated  with 
the  DA  scheme.  Error  charts  for  the  sea  surface  elevation  demonstrated  the  model’s  ability 
to  fit  the  data  within  the  error  bars  and  exposed  coastal  areas  requiring  better  coverage  by 
observations. 

Another  study  of  tidal  currents  was  made  by  Kurapov  et  al.  [24,  25]  who  retrieved  har¬ 
monic  tidal  constituents  from  the  HER  observations  off  the  Oregon  coast.  The  3d  baroclinic 
spectral  tidal  model  was  coupled  with  the  representer-based  4dVar  method  [8].  Computa¬ 
tions  with  synthetic  data  (tidal  velocity  ellipses)  showed  that  HFR  observations  at  the  sur¬ 
face  can  also  be  used  to  map  tidal  flow  at  depth.  For  the  M2  tidal  constituent,  information 
from  the  surface  was  projected  in  space-time  along  the  wave  characteristics,  thus  providing 
a  uniquely  detailed  picture  of  the  temporal  and  spatial  variability  of  internal  tide  on  the 
central  Oregon  shelf. 

In  a  more  general  study,  Iloteit  et  al.  [19]  combined  a  high  resolution  configuration 
of  the  Massachusetts  Institute  of  Technology  general  circulation  model  and  its  adjoint  to 
obtain  a  dynamically  consistent  interpolation  pattern  that  matches  HFR  data  collected  off 
the  San  Diego  coast  in  California.  The  DA  scheme  adjusted  initial  conditions,  boundary 
conditions,  and  atmospheric  forcing  fields  to  match  a  model  solution  to  observations.  The 
dynamically  constrained  interpolation  experiments  have  shown  that  the  model  solutions  can 
be  tuned  to  match  HFR  observations  within  the  error  bars  over  the  period  of  10  days,  largely 
through  the  adjustment  of  the  wind  stress  controls.  Similar  results  were  obtained  in  [38,  6] 
who  identified  wind  uncertainties  as  the  major  source  of  errors  in  coastal  models.  Thus, 
continuous  assimilation  of  HFR  data  provides  an  efficient  and  dynamically  consistent  tool 
to  correct  these  errors  and  significantly  increases  the  24-hour  forecast  skill  in  the  modeling 
of  sea  surface  currents. 

4.  Reconstruction  of  the  Missing  Data 

As  it  has  been  already  mentioned,  the  unique  feature  of  HFR  observations  is  their  nearly 
continuous  temporal  and  spatial  coverage.  However,  the  back-scattered  HFR  signals  are 
subject  to  distortions  of  artificial  and  natural  origin  which  cause  numerous  gaps  in  spatial 
coverage.  These  gaps  may  strongly  degrade  performance  of  the  interpolation  algorithms. 

The  impact  of  the  missing  radial  velocity  data  can  be  diminished  by  taking  into  account 
temporal  statistics  of  the  HFR  observations  which  allow  to  obtain  spatial  covariances  with  a 
reasonable  degree  of  accuracy.  Therefore,  an  interpolation  algorithm  in  the  space  of  radial 
velocities  can  be  thought  out  which  fills  the  gaps  using  a  truncated  EOF  decomposition  of 
the  data  covariance  matrix.  This  technique  is  widely  used  in  processing  the  SST  imagery, 
and  can  also  be  readily  applied  to  HFR  observations  [53].  In  addition,  spatial  correlations 
between  the  radial  velocities  can  be  used  to  estimate  observational  noise,  assess  its  variance, 
and,  therefore  quantify  the  cost  function  weights  of  the  2dVar  interpolation  algorithm  (4). 
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4.1.  The  Method 

The  spectral  decomposition  provides  the  following  representation  of  the  covariance  matrix: 

C  =  UAUJi  (8) 


where  U  is  a  rectangular  matrix  whose  columns  are  the  eigenvectors  ek  (empirical  orthog¬ 
onal  functions,  EOFs)  of  C  corresponding  to  the  eigenvalues  X*,  and  A  =  diag{X*}.  The 
eigenvalues  quantify  time  variation  of  the  spatial  patterns  in  the  radial  velocity  distributions 
described  by  the  corresponding  EOFs. 

Representation  (8)  can  be  employed  to  estimate  the  noise  level  using  the  cross- 
validation  (CV)  technique  (e.g.,  [7]).  The  technique  provides  a  certain  number  of  EOFs 
(modes)  Kr  that  are  well-resolved  by  the  HFR  observations.  The  rest  of  the  modes  ek,  k>  Kr 
are  attributed  to  noise,  whose  spatial  variability  cannot  be  determined  with  statistical  confi¬ 
dence  from  the  data.  Technically,  Kr  is  computed  as  the  number  of  modes  which  minimize 
the  interpolation  error  at  the  randomly  chosen  set  coc  of  CV  points.  These  points  are  tem¬ 
porarily  removed  from  observations  and  constitute  a  small  portion  of  the  data  set  in  order 
to  minimize  their  impact  on  the  result  of  covariance  estimate. 

To  fill  a  gap  containing  points  x,-  in  a  subdomain  (0  C  ft,  the  radial  velocities  v(x*) 
observed  outside  the  gap  (x,€  ft\co)  are  expanded  in  Kr  “resolved"’  eigenfunctions**  of  C: 


find  a*  : 
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and  the  expansion  coefficients  ^  are  used  to  obtain  radial  velocities  within  the  gap: 


Kr 

v(x,  e  (o)  =  e  (0)  (10) 

k  =  l 

After  this,  the  EOF  expansion  is  iteratively  improved:  a  set  of  EOFs  {ekm  }  on  the 
/nth  iteration  is  computed  using  the  covariance  estimate  C^m\  emerging  from  the  data  set 
whose  gaps  are  already  filled  with  the  help  of  the  previous  EOFs  {**„,_!)},  then  these  new 
EOFs  {**m>  }  are  employed  to  fill  the  gaps  again.  The  process  terminates  when  the  relative 
reduction  of  the  mean  interpolation  error 
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computed  over  the  CV  set  C0c  becomes  smaller  than  the  machine  precision.  Computations 
described  by  (9-11)  are  conducted  for  several  values  of  Kr  to  find  the  optimal  one  that 
minimizes  £2. 

With  the  optimal  cutoff  number  of  modes  Kr ,  the  covariance  matrix  C  can  be  decom¬ 
posed  into  the  well-resolved  Cr  and  unresolved  (noisy)  Cn  constituents: 


C  =  Cr  +  Cn  =  UrArUj  +  UnA„UTn 


(12) 
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where  \Jr  is  the  Kr  x  K  matrix,  whose  columns  are  the  first  (well-resolved)  eigenvectors, 
Ar  =  diag{X*},  k  =  l,...,tfr;  the  eigenvectors  in  the  columns  of  the  (K  —  Kr)  xK  matrix  U„ 
are  attributed  to  noise,  and  A„  =  diag{ X*  } ,  k  =  Kr  +  1 , . . . ,  K. 

The  noise  level  v  is  estimated  as 


v  = 


TrA„' 

i  ■ 

TrA 

K  K 


L  vl> 

l  =  Kr+l  /=  1 


i 


(13) 


whereas  observation  error  variances  a* ,  k  =  1 , K  (Eq.  (4))  are  represented  by  the  diago¬ 
nal  elements  of  the  matrix  C„  =  U„A„Uj.  The  diagonal  elements  of  C  can  also  be  utilized 
to  estimate  the  variances  and  c2d  of  the  respective  fields  Av,  Acurlv  and  Adivv  which 
enter  the  equation  (4)  in  the  form  of  the  inverse  variances  Wu>  Wc  and  Wd. 

In  the  next  section,  the  benefits  of  the  GF  technique  are  assessed  in  a  series  of  experi¬ 
ments  with  simulated  and  real  data. 
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Figure  8.  Setting  of  the  gap  simulation  experiments:  Gray  shading  shows  data  acquisition 
areas  of  the  radars  in  the  presence  of  a  simulated  ship  (case  b)  moving  across  the  Monter¬ 
rey  Bay.  Area  within  the  rectangle  shows  the  boundary  of  the  data-free  region  in  case  c. 
Numbers  enumerate  radars  switched  off  in  case  d. 


4.2,  Testing 

4.2.1.  Experiments  with  Simulated  Data 

Experiments  with  simulated  data  were  configured  to  mimic  real  observations  conducted 
in  the  Monterrey  Bay  in  summer  of  2003  [48]:  The  “true”  currents  u!  were  extracted 
from  the  12.5-day  run  of  the  NCOM  model  [3]  forced  by  COAMPS  [18]  winds.  Sur¬ 
face  currents  were  sampled  every  hour  along  the  beams  of  three  radars  which  probed  the 
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radial  components  v*k  of  the  model  surface  velocity  at  386,  407  and  349  points  respec¬ 
tively.  Therefore  the  dimension  of  the  data  space  was  K=  1 142.  The  total  number  of  the 
grid  points  where  velocity  vectors  were  reconstructed  was  M- 560,  so  the  number  of  un¬ 
knowns  2A/=1 120  was  approximately  equal  to  the  number  of  observations.  Overall  there 
are  N  =  Kn  =  1 142  •  301  =  343,742  observation  points,  where  n  is  the  number  of  hourly 
time  steps  in  the  12.5  day  time  window.  Observations  were  modeled  by  eq.  (5) 

Three  noise  level  values  v  (0,  0.1 ,  and  0.3)  were  tested  within  each  of  six  major  series 
of  experiments.  Each  series  was  characterized  by  specific  structure  of  the  artificial  gaps 
introduced  into  the  simulated  data  set  to  assess  the  benefits  of  the  GF  technique.  These 
simulated  data  sets  were  the  following: 

0)  without  the  gaps 

a)  with  1  -point  gaps  randomly  distributed  along  the  beams  (data  loss  y  =  1 3.5%) 

b)  with  gaps,  generated  by  obstacles,  moving  across  the  domain  (Fig.  8):  Each  obstacle 
(ship)  spoils  data  along  three  beams,  whose  intersection  point  coincides  with  the  ship’s 
position.  Back-and  -forth  motion  of  three  ships,  which  effectively  removed  6.9%  of  the 
data  points  from  observations  was  simulated. 

c)  with  the  gap  created  by  discarding  all  observation  points  in  the  rectangular  region 
(Fig.  8)  for  1  day.  This  gap  removed  28%  of  the  data  on  August  10-1 1  and  approximately 
y  =2%  of  the  data  overall. 

d)  with  gaps  generated  by  switching  off  for  1 2  hours  radars  2,3  on  August  4,  radar  1  on 
August  8  and  radars  1,3  on  August  12  (Figure  8,  T=6.3%). 

e)  with  all  the  above  mentioned  gaps  superimposed  (7^=28. 2%) 

The  quality  of  interpolation  was  monitored  by  the  same  three  parameters  eu,ec  and  ej 
as  described  by  eq.  (6-7),  The  CV  set  coc  was  specified  by  randomly  removing  10-13  points 
on  each  time  layer  with  the  total  number  of  CV  points  /Vcv=3,490  (approximately  1  %). 


Figure  9.  Interpolation  error  e  as  a  function  of  the  number  of  modes  m. 

The  cutoff  number  Kr  was  determined  for  the  noise  levels  (v=0,  0.1  and  0.3)  by  mini¬ 
mizing  the  interpolation  error  (1 1),  and  was  found  to  be  independent  on  random  variations 
in  the  space-time  structure  of  the  CV  set.  Figure  9  shows  calculations  of  Kr  for  the  experi¬ 
ments  with  v=0. 1  and  0.3:  for  observations  specified  by  eq.  (5)  the  S/N  separation  number 
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is  38  for  v=0.1  and  20  for  v=0.3.  The  corresponding  estimates  of  the  noise  level  (Eq.  (13)) 
are  0.093  and  0.29  in  very  good  agreement  with  the  true  values.  For  the  case  of  perfect 
observations  (v  =  0)  Kr  appeared  to  be  close  to  K  as  the  dependence  (11)  t(Kr)  flattened 
out  at  large  Kr  and  did  not  show  any  distinct  minimum. 

Table  3  shows  the  improvement  in  performance  of  the  2dVar  algorithm  when  the  pre¬ 
liminary  GF  is  applied. 


Table  3.  Dependence  of  the  interpolated  field  parameters  evy  ecy  and  ej  on  the  structure 
of  the  gaps  in  HFR  observations  for  v  =  0.3.  The  results  of  standard  2dVar  (without 
GF)  and  2d  Var  with  the  gaps  filled  are  shown  respectively  in  the  left  and  right  columns 
of  the  table  cells 


case 

Y 

e 

V 

Cd 

e 

0 

0.0% 

— 

0.251 

— 

0.652 

— 

0.537 

a 

13.5% 

0.260 

0.252 

0.665 

0.659 

0.551 

0.542 

b 

6.9% 

0.256 

0.251 

0.662 

0.654 

0.545 

0.540 

c 

2.0% 

0.254 

0.251 

0.660 

0.655 

0.545 

0.539 

d 

6.3% 

0.280 

0.258 

0.677 

0.661 

0.564 

0.542 

abed 

27.9% 

0.311 

0.274 

0.703 

0.679 

0.589 

0.558 

In  the  case  a  (randomly  distributed  l -point  gaps)  the  improvement  is  significantly  lower 
than  the  percentage  of  data  loss  (13.5%),  primarily  because  filling  random  1 -point  gaps 
affects  information  content  on  the  grid  scale  which  is  poorly  resolved  anyway.  Besides, 
data  absence  is  largely  compensated  by  observations  in  the  points  located  in  the  immediate 
vicinity  of  the  1 -point  gaps  at  distances  often  smaller  than  the  grid  step.  These  neighboring 
points  compensate  missing  data  and  provide  the  2dVar  interpolation  with  enough  informa¬ 
tion  on  the  larger-scale  variability.  Also  note  that  the  surface  velocity  field  is  recovered  in 
most  cases  with  a  better  accuracy  ev  than  the  noise  level  v  =  0.3. 

Moving  ships  (case  b)  spoil  6.9%  of  the  entire  set  of  343,742  data  points.  Relative 
improvement  of  the  2dVar  interpolation  (line  3  in  Table  3)  is  somewhat  smaller  than  for 
the  case  of  completely  random  gaps:  Velocity  field  is  better  by  1.1%  whereas  vorticity 
and  divergence  fields  show  0.9  and  1.0%  improvements  respectively.  Nevertheless,  the 
improvement  appears  to  be  pretty  robust  with  respect  to  this  type  of  gaps  as  well. 

Much  more  difficulties  emerge  when  a  gap  occupies  a  significant  portion  of  the  interpo¬ 
lation  domain,  as  in  case  (c)  (Fig.  8).  To  better  illustrate  the  benefits  of  the  GF  technique, 
we  placed  the  gap  at  the  location  of  an  eddy-like  structure  seen  in  the  mouth  of  the  Monter¬ 
rey  Bay  around  the  10th  of  August,  2003  (Fig.  10,  left  panel).  This  eddy  is  not  reproduced 
by  the  2d  Var  technique  alone  (Fig.  10,  right  panel),  simply  because  there  is  no  information 
on  the  eddy  in  the  velocity  field,  surrounding  the  gap.  On  the  contrary,  if  GF  is  used  prior 
to  2dVar,  a  certain  portion  of  this  eddy  emerges  from  the  EOF  statistics,  increasing  the  skill 
of  the  2dVar  algorithm  (Fig.  10,  middle  panel). 

Table  3  does  not  give  full  impression  of  the  improvement,  because  error  data  are  aver¬ 
aged  over  the  whole  observation  period  (12.5  days),  whereas  the  data  sets  in  case  (c)  differ 
only  on  the  10-1 1th  of  August.  If  averaging  is  performed  over  the  time  period  containing 
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Figure  10.  The  “true”  velocity  field  (left  panel)  and  velocity  fields  reconstructed  by  2dVar 
with  preliminary  filling  the  rectangular  gap  (middle  panel)  and  without  (right  panel). 

the  gap  (from  12PST  10.08.2003  to  12PST  1 1.08.2003)  then  the  advantage  is  obvious:  ev  is 
reduced  from  0.49  to  0.32  (33%  reduction).  Similar  error  reductions  are  also  observed  for 
the  vorticity  and  divergence  fields  (28  and  37%  respectively). 

A  common  reason  for  low  data  return  of  a  HFR  system  is  malfunction  of  one  or  more 
radars.  This  kind  of  situation  was  simulated  by  switching  off  both  northern  radars  for 
half  a  day  on  the  4th  of  August,  southern  radar  on  the  8th  and  two  southern  radars  on  the 
12th.  The  strongest  reduction  of  the  interpolation  errors  occurred  on  the  12th  of  August, 
when  the  reconstructed  (true)  currents  were  generally  perpendicular  to  the  beams  of  the 
only  operating  radar.  In  that  case  the  velocity  error  ev  reduced  64%  (from  0.84  to  0.34) 
with  66%  of  the  missing  data  being  filled.  Interpolated  velocity  distributions  show  that 
2dVar  tends  to  align  velocities  along  the  beams  of  the  only  working  radar,  producing  rather 
unrealistic  patterns.  After  filling  of  the  missing  data  from  the  southern  radars,  the  skill  of 
the  2dVar  algorithm  was  significantly  improved. 

Finally,  all  the  gaps  were  combined  together  to  obtain  a  “realistic”  HFR  record,  char¬ 
acterized  by  72%  of  the  data  return.  Figure  1 1  gives  an  overall  comparison  between  the 
methods  in  terms  of  ev t  ej  and  ec.  It  is  obvious  that  EOF-based  GF  of  the  radial  data  is 
particularly  advantageous  during  the  severe  data  loss  events  caused  by  either  malfunction 
of  a  radar  (8.8)  or  two  (4.8,  12.8);  or  by  data  loss  in  a  region,  whose  size  is  considerably 
larger  than  the  grid  step  ( 1 1 .8). 

Beyond  these  periods,  when  only  1 -point  and  ship-generated  gaps  are  present,  prelimi¬ 
nary  GF  still  has  some  (1-3%)  advantage  over  the  stand-alone  2dVar  in  terms  of  eVf  ej  and 
ec  (see  Table  3  and  Fig.  1 1). 

It  is  also  noteworthy  that  the  GF  technique  allows  to  retrieve  the  sea  surface  state  with 
the  accuracy  of  ev=0.27,  that  is  better  than  the  noise  level  v=0.3  (Fig.  1 1 )  even  in  the  case  of 
28%  loss  of  observations.  The  conventional  2dVar  technique  (ev=0.3 1 ,  blue  line  in  Fig.  1 1 ) 
demonstrates  somewhat  lower  skill,  primarily  because  of  much  poorer  performance  during 
the  heavy  data  loss  periods. 
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Figure  1 1.  Velocity  interpolation  error  ev  for  the  “perfect”  data  set  (without  gaps,  red  line), 
for  the  gappy  data  with  (black)  and  without  (blue)  preliminary  GF.  Shaded  area  denotes  the 
part  of  data  occupied  by  the  gaps.  Particularly  severe  losses  of  data  are  observed  during  sim¬ 
ulated  radar  malfunctions  (4,  8  and  12  of  August)  and  on  August  1 1th,  when  “observations” 
were  removed  from  a  large  region  shown  in  Fig.  8. 

4.2.2.  HFR  Observations  in  the  Eastern  English  Channel  and  Iroise  Sea 

To  test  the  algorithm  with  real  data,  the  GF  algorithm  was  applied  to  HFR  observations 
obtained  off  the  Opal  coast  of  France  [42]  and  in  the  Iroise  Sea  (northern  Bay  of  Biscay). 

The  Opal  coast  experiment  was  conducted  in  May-June  2003  with  two  HFRs  deployed 
on  the  Cape  Gris  Nez  (CGN)  and  Wimereux  (WMX,  Figure  12).  The  entire  35-day  record 
from  0.00  GMT  01.05.2003  to  23.40  GMT  04.06.2003  was  used  for  testing.  Surface  cur¬ 
rents  were  sampled  every  20  minutes  at  10°  azimuthal  resolution  defined  by  the  beam  width. 
The  radial  velocity  data  were  available  along  the  beams  at  1 .8  km  resolution.  Grid  cells  with 
less  than  75%  data  returns  were  excluded  from  consideration,  constraining  the  interpolation 
domain  to  the  ranges  less  than  20  km  [43]  and  the  total  number  of  quality  data  points  to 
K=203.  Overall,  the  analyzed  records  were  characterized  by  87%  of  data  return.  Approx¬ 
imately  10%  of  the  missing  radial  velocities  were  due  to  data  acquisition  problems  at  the 
Wimereux  radar  on  May  3  (3  hour  gap)  and  May  21-22  (21  hours).  The  instrumental  ac¬ 
curacy  of  the  measurements  was  5  cm/s  and  the  typical  magnitude  of  the  observed  radial 
velocities  was  Vr  =0.35  m/s. 

A  set  of  CV  points  0)c  was  generated  by  randomly  removing  8-12  radial  velocities  from 
the  data  every  20  minutes.  In  total,  24,097  (4.7%)  observations  were  removed. 

The  quality  of  inteipolation  was  estimated  as  the  mean  absolute  difference  between  the 
values  of  the  interpolated  velocity  at  the  CV  points  and  the  radial  velocities  measured  at 
these  points: 

K  =  (|v*-  Av-r*|) 

Here  index  k  enumerates  the  CV  points  and  angular  brackets  denote  averaging  over  (0C.  The 
total  number  of  gaps  in  observations  (including  the  CV  points)  was  92,156  (18%). 

Similar  to  twin-data  experiments,  the  noise  level  was  determined  by  minimizing  the  GF 
error  (11).  Dependencies  of  the  normalized  interpolation  errors  e  and  e*  on  the  number  of 
modes  were  similar  to  those  shown  in  Fig.  9  and  demonstrated  distinct  minima  at  Kr  = 
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Figure  12.  Surface  velocity  at  14.20GMT  on  May  21  2003  in  the  ERMANO  study  area. 
Contours  show  the  bathymetry  in  meters.  Gray  dots  indicate  locations  of  the  surface  veloc¬ 
ity  measurements  by  two  radars  (shown  by  circles). 


33  -  35.  The  value  Kr  =  33  was  selected  as  the  noise  cutoff  number.  The  corresponding 
observational  noise  level  computed  through  eq.  (13)  was  close  to  0.15,  or  5.1  cm/s,  in  a 
good  correspondence  with  the  above  estimate  of  the  instrumental  errors. 

Overall,  the  GF  technique  enabled  to  reduce  the  time  averaged  relative  interpolation 
error  e*v/Vr  to  0.16  (5.5  cm/s),  a  value  very  close  to  the  observational  noise  level.  On  the 
contrary,  the  mean  value  of  e*  without  preliminary  gap  filling  appeared  to  be  more  than  two 
times  higher  (0.35)  indicating  a  significant  benefit  of  combining  2dVar  interpolation  with 
the  GF  analysis. 


Figure  13.  Tidal  ellipses  computed  from  the  PCA  of  the  surface  velocities  obtained  by  the 
standard  2dVar  method  (a);  and  by  the  2dVar  method  with  GF  technique  (b,c).  The  24-hour 
PCA  averaging  is  performed  over  the  period  of  one  operating  radar  (starting  12.40GMT 
21.05)  for  (a)  and  (b)  and  over  two  24-hour  intervals  preceding  the  gap  (starting  12.40 
20.05)  and  following  immediately  after  (starting  12.00  22.05)  for  (c).  Every  second  ellipse 
is  shown. 

The  advantage  of  the  GF  technique  is  most  vividly  seen  during  the  periods  when  the 
WMX  radar  was  not  working.  As  an  example.  Fig.  12  demonstrates  a  velocity  snapshot 
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within  the  second  (21 -hour)  gap.  As  it  was  proved  by  comparison  with  the  velocity  pat¬ 
terns  obtained  for  the  identical  tidal  phase  on  May  20  and  23,  the  reconstructed  velocities 
represent  well  the  realistic  current  field.  Fig.  13  shows  tidal  ellipses  from  Principal  Com¬ 
ponent  Analysis  (PCA)  of  the  interpolated  currents  over  the  24-hour  period  from  12:00 
GMT  05.2 1 .2003  to  1 2:00  GMT  05.22.2003.  In  this  time  interval,  the  WMX  radar  was  not 
operating  for  nearly  21  hours  (from  12:40  GMT  21.05  to  8:20  GMT  22.05).  The  pattern 
in  Figure  13a  (obtained  without  preliminary  GF)  appears  to  be  completely  unrealistic  as 
the  major  axes  of  the  ellipses  tend  to  align  along  the  beam  directions  of  the  only  operating 
radar  at  Cape  Gris  Nez.  Figure  13b  is  apparently  more  close  to  reality  since  the  spatial 
distribution  of  the  ellipses  is  much  more  similar  to  the  ones  obtained  from  PCA  analyses 
of  two  24-hour  periods:  one  immediately  before  and  another  one  immediately  after  the  gap 
(Fig.  13c).  During  these  two  periods  both  radars  were  in  full  operation  with  the  average 
data  return  of  approximately  92 


Figure  14.  Right:  surface  currents  in  the  Iroise  Sea  (West  Brittany)  obtained  from  radial 
velocities  measured  by  single  radar  (not  shown,  located  5  km  east  of  SE  comer  of  the 
domain)  on  July,  23,  2008;  velocities  of  the  northern  radar  were  interpolated  using  the  GF 
technique.  Left:  surface  currents  in  the  same  area  measured  by  both  radars  two  tidal  periods 
earlier  (July  22,  2008).  50  and  100  m  isobaths  are  shown  by  gray  contours.  Northern  radar 
at  Point  Garchine  is  shown  by  black  dot. 


In  a  similar  way,  the  GF  technique  was  applied  for  processing  of  the  HFR  data  in  the 
Iroise  Sea  (section  2.2.2)  acquired  July  21  to  September  15,  2008.  During  that  period,  the 
radial  velocity  records  had  two  major  gaps.  The  first  4-hour  gap  was  caused  by  emission 
interruption  of  the  northern  radar  on  July  23  between  18:00  and  22:00GMT.  The  second, 
19-hour  long,  was  caused  by  interruption  of  the  southern  (Point  Brezellec)  radar  on  July 
28-29.  Tidal  currents  were  substantially  stronger  on  July  23  than  on  July  28  thus  inducing 
a  rotational  character  of  local  circulation  around  the  Ushant  Island  and  Sein  Archipelago. 

Figure  14  shows  snapshots  of  surface  currents  separated  by  exactly  two  tidal  cycles 
(on  July  23  and  July  22,  2008).  The  circulation  pattern  in  Fig.  14a  has  been  obtained 
from  one-radar  data  using  the  GF  algorithm  in  conjunction  with  2dVar  interpolation.  This 
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allowed  reconstruction  of  the  surface  currents  at  a  level  of  details  that  was  not  previously 
available.  An  anti-cyclonic  eddy  north  of  the  Ushant  Island,  generated  at  the  end  of  flood,  is 
fully  resolved,  general  pattern  of  circulation  during  the  current  reversal  is  well  reproduced, 
as  well  as’ the  current  intensification  in  the  Fromveur  Strait,  All  these  features  of  local 
circulation  can  be  also  found  on  the  map  derived  from  radial  velocities  recorded  by  both 
radars  one  day  before  (Fig.  14b).  Discrepancy  between  the  radial  velocities  observed  by 
the  southern  radar  and  the  interpolated  velocities  projected  on  radar  beams  remains  below 
30%.  Further  assessment  of  the  GF  algorithm  is  currently  underway  in  the  regions  with  low 
tides,  (Gulf  of  Lyons  in  the  Mediterranean  Sea  [15]). 

5.  Summary 

Monitoring  surface  currents  in  coastal  areas  with  HFRs  is  a  rapidly  developing  application 
of  the  radar  technology.  Of  particular  importance  is  the  HFRs’  capability  to  observe  circula¬ 
tion  up  to  150-200  km  offshore  with  spatial  and  temporal  resolutions  sufficient  for  accurate 
tracking  and  prediction  of  floating  particle  trajectories,  be  it  oil  spills  or  rafts  with  people 
waiting  to  be  rescued. 

In  this  chapter  we  have  reviewed  the  latest  developments  in  the  algorithms  for  retriev¬ 
ing  the  surface  velocity  field  from  the  radial  velocities.  The  radials  are  derived  by  process¬ 
ing  backscattered  HFR  signals,  a  technique  equally  important  for  increasing  spatial  resolu¬ 
tion  and  overall  quality  of  the  surface  velocity  field  (see,  for  example,  [31]  and  references 
therein). 

In  Section  2,  a  numerical  algorithm  based  on  minimization  of  a  quadratic  cost  function 
in  the  space  of  all  possible  configurations  of  the  velocity  field  was  presented.  Performance 
of  the  method  is  compared  with  the  LI  and  OMA  algorithms.  It  is  shown  that  the  2dVar 
approach  is  robust  and  capable  to  provide  a  statistically  consistent  fit  to  the  data  in  the 
wide  range  of  signal/noise  ratios.  The  comparison  demonstrated  similar  (to  LI)  or  better 
(than  OMA)  performance  of  the  2dVar  technique  under  relatively  high  signal/noise  ratios, 
especially  when  a  80%-90%  fit  to  the  velocity  field  containing  strong  localized  features 
is  required.  At  more  realistic  (less  than  3-4)  signal/noise  ratios  the  OMA  and  2dVar  have 
similar  skill  and  outperform  LI  because  of  their  better  treatment  of  the  coastal  regions  where 
beam  directions  are  close  to  parallel. 

Comparison  of  the  algorithms  has  also  shown  that  the  variational  approach  gives  more 
flexibility  in  fitting  the  data,  since  the  number  of  modes  in  use  is  close  to  the  number  of 
points  on  the  interpolation  grid.  2dVar  method  is  also  flexible  in  the  choice  of  regularization, 
because  the  desired  smoothness  and  its  spatial  variation  can  be  directly  controlled  by  a 
simple  modification  of  the  cost  function  weights.  The  method  also  appears  to  be  more 
accurate  than  LI  and  OMA  in  the  regions  with  sparse  data  coverage,  and  within  the  large 
gaps  in  observations. 

The  2dVar  algorithm  is  regularized  by  enforcing  smoothness  in  the  interpolated  patterns 
of  divergence  and  vorticity.  This  formulation  attempts  to  retrieve  spatial  distribution  of 
these  quantities  because  they  are  extremely  important  for  successful  tracking  of  surface 
contaminants  and  other  applications  (search  and  rescue,  route  optimization). 

Section  3  provides  an  overview  of  more  sophisticated  interpolation  techniques,  which 
synthesize  HFR  data  with  other  observations  and  dynamical  information  provided  by  the 
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numerical  models  of  oceanic  circulation.  This  type  of  interpolation  has  been  under  exten¬ 
sive  development  in  the  last  decade,  but  HFR  data  were  rarely  considered  even  in  regional 
applications  possibly  because  of  their  limited  availability  compared,  for  example,  to  satel¬ 
lite  observations.  Nevertheless,  it  has  been  shown  that  using  HFR  data  greatly  improves 
the  forecast  skill  of  the  DA  systems,  especially  when  increased  predictability  of  the  surface 
currents  is  required. 

An  open  question  in  HFR  data  assimilation  is  whether  the  radial  velocities  require  pre¬ 
processing  or  should  be  assimilated  directly  into  the  numerical  models.  We  assume  that 
in  many  cases  preliminary  preprocessing  could  be  beneficial.  In  particular,  kinematically 
constrained  interpolation  of  the  HFR  radials  on  the  model  velocity  grid  may  inhibit  spuri¬ 
ous  fast  waves  arising  in  the  case  of  strongly  intermittent  data  availability  which  is  often 
observed  in  the  practice  of  HFR  measurements. 

In  Section  4  we  demonstrated  the  benefits  of  more  sophisticated  preprocessing,  which 
combines  the  EOF  analysis  with  the  2dVar  interpolation  technique.  This  approach  is  able 
to  successfully  process  occasional  single-radar  coverage  events  and  improve  the  overall 
quality  of  monitoring  of  sea  surface  currents  by  the  HF  radars.  The  EOF  analysis  of  the 
radial  velocities  provides  a)  statistically  rigorous  estimation  of  the  weights  for  the  2dVar 
algorithm,  and  b)  a  set  of  spatial  patterns  (EOFs)  capable  to  fill  large  gaps  in  the  data  caused 
by  radar  malfunctions.  The  approach  takes  the  advantage  of  the  frequent  time  sampling  by 
the  HFRs  and  employs  observation  history  to  estimate  the  leading  modes  of  variability  of 
the  radial  velocities.  Numerical  experiments  with  simulated  and  real  data  have  shown  that 
preliminary  gap-filling  is  extremely  beneficial  during  occasional  periods  of  heavy  data  loss 
associated  with  radar  malfunctioning:  With  the  proposed  technique,  the  2dVar  interpolation 
errors  during  these  periods  are  typically  reduced  1.5-2  times  providing  much  more  realistic 
velocity  distributions  (Fig.  11,  14). 

Overall,  an  advanced  interpolation  method  can  be  summarized  as  a  four-step  procedure: 
EOF  analysis  of  the  radial  velocities;  estimation  of  the  noise  and  the  2dVar  cost  function 
weights;  filling  large  gaps  in  observations,  and  finally,  retrieving  of  the  velocity  vectors  from 
the  filled  data  set.  This  type  of  interpolation  eliminates  disruptions  in  HFR  observations  of 
surface  currents  caused  by  environmental  factors  and  by  discontinuities  in  HFR  operation. 

Further  development  of  the  preprocessing/interpolation  techniques  for  the  HFR  radial 
velocities  and  their  assimilation  into  the  numerical  models  is  a  necessary  prerequisite  of 
continuous  improvement  in  monitoring  near-coastal  circulation  which  is  extremely  impor¬ 
tant  in  practical  applications  and  may  help  to  solve  many  environmental  problems  caused 
by  human  activity. 


References 

[1]  Alvera-Azcarate  A.,  A.  Barth,  M.  Rixen  and  J.M.  Beckers,  2005:  Reconstruction  of 
incomplete  oceanographic  data  sets  using  empirical  orthogonal  functions:  application 
to  the  Adriatic  Sea  surface  temperature,  Ocean  Modelling ,  9,  325-346. 

[2]  Alvera-Azcarate,  A.,  A.  Barth,.  J.-M.  Beckers,  and  R.  H.  Weisber,  2007:  Multivari¬ 
ate  reconstruction  of  missing  data  in  sea  surface  temperature,  chlorophyll,  and  wind 
satellite  fields,  J.  Geophys.  Res,,  1 12,  C03008,  doi:  10.1029/2006JC003660. 


Interpolation  of  the  Radial  Velocity  Data ... 


163 


[3]  Barron,  C.N.,  A.B.  Kara,  H.E.  Hurlburt,  C.  Rowley,  and  L.F.  Smedstad  (2004),  Sea 
surface  height  predictions  from  the  global  Navy  Coastal  Ocean  Model  (NCOM)  dur¬ 
ing  1998-2001  ,J.  Atmos.  Oceanic  Technol.,  21(12),  1876-1894. 

[4]  Barth,  A.,  Alvera-Azcrate,  A.,Weisberg,  R.H.,  2008.  Assimilation  of  high-frequency 
radar  currents  in  a  nested  model  of  the  West  Florida  Shelf.  J.  Geophys.  Res.,  1 13, 
C08033,  doi.  10. 1029/2007JC004585. 

[5]  Barth,  A.,  Alvera-Azcrate,  A.,  Beckers,  J.M.,  Staneva,  J.,  Stanev,  E.,  Schulz- 
Stellenflcth,  J.,  2010:  Correcting  surface  winds  by  assimilating  High-Frequency  Radar 
surface  currents  in  the  German  Bight.  Ocean  Dynamics,  61,599-610. 

[6]  A.  Barth,  A.  Alvera-Azcrate,  K  -W.  Gurgel,  J.  Staneva,  A.  Port,  J.-M.  Beckers,  and 
E.  V.  Stanev,  2010:  Ensemble  perturbation  smoother  for  optimizing  tidal  boundary 
conditions  by  assimilation  of  High-Frequency  radar  surface  currents  -  application  to 
the  German  Bight.  Ocean  Science,  6(1),  161 178,  doi:  10.5 194/os-6-l 6 1-2010. 

[7]  Beckers,  J.M.  and  M.  Rixen,  2003:  EOF  calculations  and  data  filling  from  incomplete 
oceanographic  observations,  J.  Atm.  Oceanic  Tech.,  20(1 2),  1839- 1856. 

[8]  Bennett,  A.F.,  2002:  Inverse  modeling  of  the  ocean  and  atmosphere.  Cambridge  Uni¬ 
versity  Press,  256  pp 

[9]  Breivik,  O.,  and  O.  Ssetra,  2001:  Real  time  assimilation  of  HF  radar  currents  into  a 
coastal  ocean  model,  J.  Marine  Sys.,  28,161-182. 

[10]  Chapman,  R.  D.,  L.  K.  Shay,  H.  C.  Graber,  J.  B.  Edson,  A.  Karachintsev,C.  L.  Trump, 
and  D.  B.  Ross,  1997:  On  the  accuracy  of  the  HF  radar  surface  measurements:  Inter- 
comparisons  with  ship-based  sensors.  J.  Geophys.  Res.,  102,  18,737-18,748. 

[11]  C.  Chavanne,  1.  Janekovic,  P.  Flament,  P.-M.  Poulain,  M.  Kuzmic,  and  K.-W. 
Gurgel,  2007:  Tidal  currents  in  the  northwestern  Adriatic:  High-frequency  radar 
observations  and  numerical  model  predictions,  J.  Geophys.  Res.,  112,  C03S21, 
doi:  10. 1029/2006JC003523. 

[  1 2]  Devenon,  J-L.,  1 990:  Optimal  control  theory  applied  to  an  objective  analysis  of  a  tidal 
current  mapping  by  HF  Radar.  J.  Atmosph.  Oceanic  Tech.,  7,  269-284. 

[13]  Eremeev,  V.N.,  L.M. Ivanov  and  A.D.  Kirwan,  1992:  Reconstruction  of  oceanic  flow 
characteristics  from  quasi-Lagrangian  data:  1.  Approach  and  mathematical  methods. 
J.  Geophys.  Res.,  97(C6),  9733-9742. 

[14]  Fang,  F„  C.C.  Pain,  I.  M.  Navon,  M.D.  Piggott,  G.J.  Gorman,  P.E.  Farrell,  P.A.  Alli¬ 
son,  and  A.J.H.  Goddard,  2008.  A  POD  reduced-order  4D-var  adaptive  mesh  ocean 
modelling  approach,  7/i/.  J.  Nutn.  Methods  in  Fluids,  doi:20.2002/fld.l91 1. 

[15]  Forget,  P,  Y.  Barbin,  and  G.  Andre,  2008:  Monitoring  of  surface  ocean  circulation 
in  the  Gulf  of  Lions  (North  West  Mediterranian  Sea)  using  WERA  HF  radars,  Proc. 
1GARSS,  CD  ROM,  Boston,  USA  (7-1 1  July). 


164 


Max  Yaremchuk  and  Alexei  Sentchev 


[16]  Graber,  H.  G,  D.  R.  Thompson,  and  R.  E.  Carabde,  1996:  Ocean  surface  features 
and  currents  measured  with  synthetic  aperture  radar  interferometry  and  HF  radar,  J. 
Geophys.  Res .,  101(C1 1),  25,813-25,832. 

[17]  Hisaki,  Y.,  T.  Tokeshi,  W.  Fujie,  K.  Sato,  and  S.  Fujii,  2001 :  Surface  current  variability 
east  of  Okinawa  Island  obtained  from  remotely  sensed  and  in  situ  observational  data, 
J.  Geophys.  Res.,  106(C12),  31,057-31,073. 

[18]  Hodur,  R.  M„  J.  Pullen,  J.  Cummings,  J„  X.  Hong,  J.D.  Doyle,  RJ.  Martin, 
M.A.Rennick,  2002:  The  coupled  ocean-atmospheric  mesoscale  prediction  system 
(COAMPS),  Oceanography ,  15(1),  88-98. 

[19]  Hoteit,  1  ,  B  Comuelle,  S.Y.  Kim,  G.  Forget,  A.  Kohl  and  E.  Terrill,  2009:  Assessing 
4D-VAR  for  dynamical  mapping  of  coastal  high-frequency  radar  in  San  Diego,  Dyn. 
Atmosph.  Oceans ,  48,  175-197. 

[20]  Kaplan,  D.,  and  F.  Lekien,  2007:  Spatial  interpolation  of  surface  current  data 
based  on  open-boundary  modal  analysis.  J.  Geophys.  Res.,  112,  C12007,  doi: 
10. 1029/2006JC003984. 

[21]  Kim,S.  Y.,  E.  Terrill,  and  B.  Comuelle,  2007.  Objectively  mapping  HF  radar-derived 
surface  current  data  using  measured  and  idealized  data  covariance  matrices,  J.  Geo¬ 
phys.  Res ,  1 12,  C0602 1 ,  doi:  1 0. 1 029/2006JC003756. 

[22]  Kim,S.  Y.,  E.  Terrill,  and  B.  Comuelle,  2008.  Mapping  surface  currents  from  HF 
radar  radial  velocity  measurements  using  optimal  interpolation,  J.  Geophys.  Res ,  1 13, 
Cl 0023,  doi:  10.1 029/2007 JC004244. 

[23]  Kondrashov,  D.  and  M.  Ghil,  2006:  Spatio-temporal  filling  of  missing  points  in  geo¬ 
physical  data  sets,  Nonlin.  Processes  Geophys.,  13,151-159. 

[24]  Kurapov,  A.  L.,  G.D.  Egbert,  J.S.  Allen,  and  R.N.  Miller,  2009:  Representer-based 
analyses  in  the  coastal  upwelling  system,  Dyn.  Atmosph.  Oceans,  48,  198-218. 

[25]  Kurapov,  A.L.,  G.D.  Egbert,  J.  S.  Allen,  R.  N.  Miller,  S.Y.  Erofeeva,  P  M.  Kosro, 
2003:  The  M2  Internal  Tide  off  Oregon:  Inferences  from  Data  Assimilation.  J.  Phys. 
Oceanogr.,  33,  17331757. 

[26]  Lekien,  F.,  and  C.  Coulliette,  2004:  Open-boundary  modal  analysis:  Interpolation,  ex¬ 
trapolation  and  filtering.  J.  Geophys.  Res.,  109,C12004,doi:10.1029/2004JC002323. 

[27]  Le  Provost  C.  and  A.  Poncet,  1978:  Finite  element  method  for  spectral  modelling  of 
tides,  bit.  J.  Num.  Meth.  in  Engineer.,  12,  853-871. 

[28]  Lewis,  J.K.,  1  Shulman,  and  A.F.  Blumberg,  1998:  Assimilation  of  the  Doppler  radar 
current  data  into  numerical  ocean  models,  Cont.  Shelf  Res.,  18,  541-559. 

[29]  Li,  Z.,  Y.  Chao,  J.C.  McWilliams,  and  K.  Ide,  2009:  A  three-dimensional  variational 
data  assimilation  scheme  for  the  Regional  Ocean  Modeling  System,  J.  Atmosph. 
Oceanic  Tech.,  25,2074-2090. 


Interpolation  of  the  Radial  Velocity  Data  ... 


165 


[30]  Lipa,  B.J.,  and  D.E.  Barrick,  1983:  Least-squares  methods  for  the  extraction  of  surface 
currents  from  CODAR  crossed-loop  data:  Application  at  ARSLOE  IEEE  J.  Ocean. 
Eng .,  8(4),  226253. 

[31]  Lipa,  B.,  Nyden,  B.,  Ullman,  D.S.,  and  Terrill,  E.,  2006:  SeaSonde  radial  velocities: 
derivation  and  internal  consistency.  IEEEJ.  Oceanic  Eng.,  31(4),  850-861. 

[32]  Lipphardt,  B.  L.,  A.  D.  Kirwan,  C.  E.  Grosch,  J.  K.  Lewis,  and  J.  D.  Paduan,  2000. 
Blending  HF  radar  and  model  velocities  in  Monterey  Bay  through  normal  mode  anal¬ 
ysis.  J.  Geophys.  Res.,  105,  C2,  3,425-3,450. 

[33]  Manette,  V.,  N.  Thomas,  V.  Cochin,  Y.  Guichoux,  and  F.  Ardhuin,  2006:  Premiers 
r^sultats  de  lexpdrience  SURLITOP  (Surveillance  Littorale  Op^rationnelle),  Naviga¬ 
tion ,  54,4557. 

[34]  McIntosh,  P.  C.,  1990.  Oceanographic  data  interpolation:  Objective  analysis  and 
splines.  J.  Geophys.  Res.,  95(C8),  13,529-13,541. 

[35]  Oke,  P.R.,  Allen,  J.S.  Miller,  R.N.,  and  Egbert  G.D.  2002:  A  modelling  study  of 
the  three-dimensional  continental  shelf  circulation  off  the  Oregon.  Part  I.  Dynamical 
analysis.  J.  Phys.  Oceanogr.,  32,  1383-1403.  J.  Geophys.  Res.  109,  C07S09. 

[36]  Oke,  PR.,  Allen,  J.S.,  Miller,  R.N.,  Egbert,  G.D.,  Kosro,  P.M.,  2002:  Assimilation  of 
surface  velocity  data  into  a  primitive  equation  coastal  ocean  model.  J.  Geophys.  Res., 
doi :  1 0. 1 029/2000JC005 1 1 . 

[37]  Oke,  P,  Allen,  J,  Miller,  RM  Egbert,  G.,  Kosro,  P,  2010.  Assimilation  of  surface 
velocity  data  into  a  primitive  equation  coastal  ocean  model.  Ocean  Dynamics,  60, 
1497-1537. 

[38]  Paduan,  J.D.,  Shulman,  1.,  2004:  HF  radar  data  assimilation  in  the  Monterey  Bay  area. 
J.  Geophys.  Res,  109,  C07S09. 

[39]  Park,  K.,  Z.  Li,  J.  Farrara  and  Yi  Chao,  2007:  Assimilation  of  high-frequency 
radar  measurements  for  coastal  currents  using  ROMS,  Proc.  SPIE,  vol.  6685, 
doi:  10. 1 1 17/12.759504, 55-68. 

[40]  Prandle,  D.,  1993:  Year-long  measurements  of  flow  through  the  Dover  Strait  by  HF 
radar,  Oceanologica  Acta,  16,457-468. 

[41]  Sentchev,  A.  and  M.  Yaremchuk,  1999:  Tidal  motions  in  the  Dover  Straits  as  a  varia¬ 
tional  inverse  of  the  sea  level  and  surface  velocity  data,  Cont.  Shelf  Res.,  19,  1905- 
1932. 

[42]  Sentchev,  A.,  M.  Yaremchuk,  and  F.  Lyard  2006:  Residual  circulation  in  the  English 
Channel  as  a  dynamically  consistent  synthesis  of  shore-based  observations  of  sea  level 
and  currents,  Cont.  Shelf  Res.,  26,  1884-1904. 

[43]  Sentchev,  A.  and  M.  Yaremchuk,  2007:  VHF  radar  observations  of  surface  currents 
off  the  northern  Opal  coast  in  the  eastern  English  Channel,  Cont.  Shelf  Res.,  27, 
2449-2464. 


166 


Max  Yaremchuk  and  Alexei  Sentchev 


[44]  Sentchev,  A.,  P.  Forget,  Y.  Barbin,  and  M.  Yaremchuk,  2012:  Surface  circulation  in 
the  Iroise  Sea  (Western  Brittany)  derived  from  high  resolution  current  mapping  by  HF 
radars.  J.  Marine  Systems  (in  press). 

[45]  Shay,  L.  K.,  J.  Martinez-Pedala,  T.  M.  Cook,  and  B.  K.  Haus,  2007:  High-frequency 
radar  mapping  of  surface  currents  using  WERA.  J.  Atm .  Oceanic  Tech .,  24, 484-503. 

[46]  Shulman,  I.,Wu,  C.R.,  Lewis,  J.K.,  Paduan,  J.D.,  Rosenfeld,  L.K.,  Kindle,  J.C.,  Ramp, 
S.R.,  Collins,  C.A.,  2002.  High  resolution  modeling  and  data  assimilation  in  the  Mon¬ 
terey  Bay  area.  Cont.  Shelf  Res.  22(8),  1  1291151. 

[47]  Shulman,  I.,  Kindle,  J.,  Martin,  P,  deRada,  S.,  Doyle,  J.,  Penta,  B.,  Ander¬ 
son,  S.,  Chavez,  F.,  Paduan,  J.,  Ramp,  S.,  2007.  Modeling  of  upwelling/relaxation 
events  with  the  Navy  Coastal  Ocean  Model.  J.  Geophys.  Res.  112,  C06023, 
doi:  1 0. 1029/2006JC003946. 

[48]  Shulman,  I.  and  J.  Paduan,  2009.  Assimilation  of  HF  radar-derived  radials  and  total 
currents  in  the  Monterrey  Bay  area,  Deep-Sea  Res .  //,  56,  (3-5),  149-160. 

[49]  Tipett,  M.K.,  J.  Anderson,  C.  Bishop,  T.  M.  Hamill  and  J.  S.  Whitaker,  2003:  Ensem¬ 
ble  square  root  filters.  Mon.  Wea.  Review ,  131,  1485-1490. 

[50]  Wilkin,  J.L.,  H.G.  Arango,  D.B.  Haidvogel,  C.S.  Lichtenwalner,  S.M.  Glenn,  and  H.S. 
Hedstrm,  2005:  A  Regional  Ocean  Modeling  System  for  the  long-term  ecosystem 
observatory.  J.  Geophys.  Res,  110,  C06S91. 

[51]  Yaremchuk,  M.,  and  A.  Sentchev,  2009:  Mapping  radar-derived  sea  surface  currents 
with  a  variational  method,  Cont.  Shelf  Res.,  29,1711-1722. 

[52]  Yaremchuk,  M.,  D.  Nechaev,  and  G.  Panteleev,  2009:  A  method  of  successive  correc¬ 
tions  of  the  control  subspace  in  the  reduced-order  variational  data  assimilation,  Mon. 
Wea.  Rev.,  137(9),  2966-2978. 

[53]  Yaremchuk,  M.,  ,  and  A.  Sentchev,  2011:  A  combined  EOF/variational  approach  for 
mapping  radar-derived  sea  surface  currents,  Cont.  Shelf  Res.,  31,  758-768. 


